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2.  Background 


Underground  facilities  (UGF)  may  be  detected  and  monitored  by  either  active  or  passive 
methods.  Active  methods  introduce  some  form  of  energy  into  the  ground,  typically  mechanical 
or  electromagnetic  waves.  UGF  will  be  manifested  in  this  data  since  they  are  air-filled  and,  as 
such,  present  a  substantial  change  in  mechanical  or  electrical  properties.  A  3-D  image  can  be 
reconstructed  from  this  data  that  shows  the  location  and  extent  of  a  UGF  (Witten,  Won,  and 
Norton,  1997  a  and  b;  Witten,  Schatzberg,  and  Devaney,  1996;  Witten,  1991;  Witten,  Levy, 
Ursic,  and  White,  1995). 


In  spite  of  the  success  of  these  methods,  the  active  methods  do  have  restrictions  that  limit  then- 
applicability.  First,  the  deployment  of  active  sources,  in  many  situations,  may  not  be  possible. 
Furthermore,  the  methods  require  that  sources  and  receivers  be  distributed  over  some  regular 
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grid  system.  Such  requirements  may  be  impossible  to  satisfy,  particularly  if  sensors  must  be 
dropped  from  aircraft.  There  will  always  be  a  loss  of  signal  with  distance  from  the  source  so 
that  active  systems  employing  both  sources  and  receivers  on  the  ground  surface  must  travel  a 
distance  that  is,  at  a  minimum,  twice  the  target  depth.  This  attenuation  will  always  limit  the 
depth  to  which  active  systems  can  be  applied. 

Passive  methods  detect  and  characterize  UGF  from  sources  generated  by  the  facility  itself. 
When  a  UGF  is  quiescent,  it  cannot  be  detected  by  a  passive  method.  However,  when  active,  it 
can  be  detected  by  remotely  monitoring  emissions.  Applying  specific  signal  processing 
techniques  allows  the  source  to  be  localized  and  characterized  by  its  emitted  spectrum.  These 
techniques  do  not  require  dense  or  uniform  sensor  spacing.  Since  the  energy  must  only  travel 
one  way,  upward  from  the  source  to  the  detectors,  it  is  far  less  limited  than  active  methods. 

In  time-exposure  photography,  low  intensity  light  is  integrated  over  time  to  enhance  image 
brightness  in  poor  lighting  conditions.  A  similar  concept  is  possible  in  acoustics,  known  as 
Time-Exposure  Acoustics  (TEA:  Norton  and  Won,  2000:  attached  here  as  Appendix  B), 
wherein  acoustic  images  are  reconstructed  from  passively  acquired  seismic  noise,  such  as  that 
generated  inside  a  UGF.  Noise  sources  include  motor-driven  machinery  (generators, 
ventilation  fens,  etc.)  and  moving  vehicles  in  the  UGF  or  through  tunnels  that  may  be  a  part  of 
its  infrastructure.  Figure  1  illustrates  the  creation  of  a  wave  front  impinging  on  a  surface 
geophone  array  due  to  a  moving  vehicle  inside  a  tunnel.  The  data  are  recorded  by  a  geophone 
array  on  the  surface  and  processed  on  site  in  real  time.  The  imaging  of  subsurface  sources 
using  passive  seismic  data  is  potentially  attractive  for  several  reasons: 


Figure  1:  Wavefront  impinging  on  a  geophone  array  generated  by  subsurface  vehicle. 
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•  The  TEA  allows  a  long  time  exposure  that  accumulates  faint  sources,  progressively 
improves  the  resolution,  and  displays  the  current  image  in  real  time. 

•  The  signal  is  integrated  over  time  whenever  the  source  is  "on,"  creating  an  image  that 
shows  collective  and  historic  activities.  In  this  way,  weak  sources  should  eventually 
reveal  themselves  after  a  sufficient  exposure.  For  a  moving  source,  such  as  a  vehicle 
moving  in  a  tunnel,  the  time  exposure  will  trace  the  path  of  the  vehicle. 

•  Passive  seismic  data  can  be  acquired  more  discretely  than  active-source  data,  particularly 
in  a  hostile  environment  or  when  the  use  of  active  sources  is  tactically  impractical. 

•  Since  the  image  is  constructed  in  real  time,  there  is  no  need  of  bulky  data  storage. 


3.  Simulation  Studies 

Appendix  A  describes  the  TEA  algorithm  originally  developed  by  Geophex  in  2000  (Norton 
and  Won,  2000).  We  first  apply  the  TEA  algorithm  to  synthetic  data  in  order  to  assess  its 
efficacy  and  to  identify  conditions  under  which  its  performance  could  be  degraded.  The  model 
includes  an  underground  acoustic  source  and  an  array  of  receivers  distributed  on  the  ground 
surface  directly  above  the  source.  The  simulated  sources  have  included  impulsive  points, 
incrementally  moved  impulsive  points,  and  continuous  random  and  periodic  point  sources. 

In  the  first  test,  an  impulsive  point  source  is  located  somewhere  in  the  subsurface  and  synthetic 
data  is  generated  for  a  linear  array  of  64  receivers  directly  above  this  point  (Fig.  2a).  An 
arbitrary  time  shift  is  introduced  (Fig.  2b)  and  the  TEA  algorithm  is  applied.  Figure  2c  presents 
a  2-D  display  of  the  likelihood  of  source  location  in  false  color  with  red  corresponding  to  the 
highest  computed  values.  The  dark  red  node  in  this  image  occurs  at  the  precise  source  location. 


(a)  (b)  (c) 


Figure  2:  (a)  Synthetic  data  for  an  impulsive  point  source,  (b)  that  has  been  shifted, 
and  (c)  the  image  resulting  from  the  application  of  the  TEA  algorithm. 

Figure  3a  displays  a  random  time-series  and  Fig.  3b  shows  the  synthetic  data  acquired  for  a 
subsurface  point  source  having  this  output.  The  TEA  image  is  given  in  Fig.  3  c  displaying  a 
result  comparable  to  the  impulsive  point  source  (Fig.  2c). 
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Figure  3:  A  random  time  series  (a),  the  synthetic  data  resulting  when  these  temporal 
variations  are  applied  to  a  point  source  (b),  and  (c)  the  TEA  image. 


The  final  2D  example  is  for  multiple  impulsive  point  sources.  In  this  simulation,  five  point 
sources  are  distributed  over  a  uniform  depth  in  the  subsurface  and  each  source  is  sequentially 
fired  such  that  there  is  a  time  delay  that  increases  with  horizontal  source  location.  The 
synthetic  data  is  provided  in  Fig.  4a  and  the  TEA  image  in  Fig.  4b.  If  this  reconstruction  were 
perfected,  five  distinct  highs  should  appear  at  uniform  horizontal  intervals  and  at  the  same 
depth.  These  features  do  appear  in  the  image  along  with  several  weaker  spurious  responses. 
These  artifacts  are  associated  with  locations  in  the  data  where  there  are  high  signal  strengths 
associated  with  points  of  intersection  of  individual  hyperbolas. 


(a)  (b) 

Figure  4:  Five  sequentially  fired  impulsive  point  sources  (a)  and  (b)  the  TEA  image. 


4.  Confirmatory  Field  Tests 

To  confirm  the  TEA  algorithm,  we  acquired  field  data  over  a  small  cave  in  limestone  in 
southeastern  Oklahoma.  The  data  were  acquired  over  a  30  by  35  ft  planar  and  slightly  sloping 
region  over  a  short  cave  accessed  through  an  open  sinkhole  (Fig.  5a).  An  array  of 48, 40Hz 
geophones  was  distributed  as  an  8  by  6  array  as  shown  in  Fig.  5b. 

An  acoustic  signal  was  created  within  the  cave  by  gently  tapping  on  the  limestone  surface  of 
the  cave.  Source  locations  were  randomly  selected  within  the  cave  over  about  25  ft  of  its 
accessible  length.  About  ten  blows  were  applied  to  each  source  position. 
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Figure  5:  (a)  The  sinkhole  cave  access  in  the  foreground  and  orange  pin  flags  at  geophone 
positions  in  the  background;  (b)  general  geometry  of  the  geophone  layout  relative  to  the  cave. 


5.  TEA  Applied  to  Field  Data 

Figure  6a  presents  data  acquired  from  one  particular  hammer  blow  within  the  cave.  The 
vertical  axis  is  the  time  sample  number  and  the  horizontal  coordinate  is  the  receiver  index  as 
identified  in  Fig.  5b.  A  sequence  of  hyperbolas  is  evident  in  the  data  and  these  are  from  the 
lines  of  geophones  closest  to  a  point  on  the  ground  surface  directly  above  the  source.  The 
result  of  the  application  of  the  TEA  algorithm  is  given  in  Fig.  6b  where  the  maximum 
computed  source  likelihood  is  displayed  as  a  surface  in  3-D  enclosed  in  red.  The  location  of 
the  impulsive  point  source  is  well  resolved. 


(a)  (b) 


Figure  6:  The  data  acquired  above  an  impulsive  point  source  located  within  a  cave 
(a)  and  (b)  the  source  location  as  provided  by  the  TEA  algorithm. 
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Figure  7  presents  acquired  data  similar  to  that  shown  in  Fig.  6  but  for  a  different  source 
location  and  here  the  character  of  the  data  is  distinctly  different.  Most  significant  is  the 
“ringing”  that  is  absent  in  Fig.  6.  This  effect  is  believed  to  be  caused  by  the  excitation  of 
resonant  modes  of  the  cave.  The  application  of  the  TEA  algorithm  to  an  early  time  portion  of 
this  data  set  is  provided  in  Fig.  8a  and  the  entire  cave  is  revealed  from  a  single  source  point. 
Using  the  full  data  series,  the  result  of  the  TEA  algorithm  (Fig.  8b)  clearly  shows  that  the 
resonance  apparent  in  Fig.  8  is  a  result  of  the  excitation  of  resonant  modes  of  both  the  known 
cave  and  deeper  unknown  void. 
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Figure  7:  Data  acquired  from  a  source  point  located  deep  within  the  cave. 
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Figure  8:  TEA  images  derived  from  the  data  given  in  Fig.  8  for  (a)  a  time- windowed 
subset  of  the  data  and  (b)  the  full  data  set.  Note  that  the  image  indicates  an  additional 
bigger  cave  below  the  small  known  one. 
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6.  Moving  Source  Field  Studies 

To  demonstrate  the  TEA  algorithm  for  a  moving  source,  data  were  acquired  on  the  Geophex 
property  in  Raleigh,  NC.  Specifically,  two  perpendicular  lines  of  24  geophones  each  were  used 
and  a  moving  source  was  created  by  striking  the  ground  with  a  hammer  at  various  locations. 
Figure  9  presents  the  source  locations  (circles)  and  the  two  geophone  lines.  Along  each  line, 
the  geophone  spacing  was  5  ft.  and  temporal  samples  were  acquired  at  2  ms  intervals. 
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Figure  9.  Measurement  geometry  used  for  the  moving  source  study. 

Figure  10a  shows  the  estimated  location  (red)  for  a  single  source  location  and  Fig.  10b  is  the 
results  of  the  time  exposure  over  all  data.  Here,  the  T  pattern  of  the  source  locations  is  clearly 
reconstructed. 


Figure  10.  TEA  image  of  (a)  a  single  location  and  (b)  the  time 
exposure  image  of  the  moving  source. 


7.  Phase  I  Summary  Conclusions 

A  Time  Exposure  Acoustics  algorithm  has  been  developed  for  locating  underground  noise 
sources  by  passive  monitoring  on  the  ground  surface.  This  algorithm  has  been  tested  on  both 
synthetic  data  and  real  data  with  equally  good  results.  Of  particular  importance  is  the  discovery 
that  an  entire  underground  void  system  can  be  imaged  based  on  noise  from  a  single  location. 
This  is  a  result  of  the  excitation  of  resonant  modes  within  this  subsurface  structure. 
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The  computational  performance  of  the  algorithm  is  such  that  it  is  possible  for  it  to  be  used  for 
real  time  acquisition  and  processing  for  detecting  and  monitoring  underground  facilities. 


8.  Plans  for  Phase  n 

Phase  n,  if  granted,  will  included  the  following  works  in  theory  and  field  tests: 

•  Refinement  of  the  TEA  algorithms  -  effect  of  varying  acoustic  speed  in  earth,  mode 
conversions,  and  surface  waves  (ground  rolls);  increased  image  resolution. 

•  Planned  field  tests  for  moving  vehicles  and  other  targets: 

1 .  Open  highways  -  visible  ground  truths 

2.  Highway  tunnels  -  example:  1-40  across  the  Appalachian  mountain  in  the  western  NC 

3.  Underground  mines  -  example:  coalmines  in  the  mountain;  excavators,  loaders, 

conveyer  belts,  ore  bucket  lifts,  and  occasional  blasts. 

4.  Known  active  UGFs  -  to  be  coordinated  with  AFOSR 

•  Design  study  for  a  dedicated  TEA  hardware  and  software  -  DSP  level  parallel  processing 
for  realtime  deployment 
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Appendix  A:  TEA  Imaging  Algorithm 


1.  Formulation  of  the  Problem 

Let  u(rn,t)  denote  the  signal  recorded  at  the  point  r„  on  the  surface  (i.e.,  r„  is  the 
position  of  the  n-th  geophone)  and  let  c  be  the  velocity  of  sound  (assumed 
constant).  See  figure  1. 


The  image,  0(r),  is  computed  by  backpropagating  the  received  signals,  u(rn,  t), 
back  to  the  source  point,  r,  and  summing  over  the  geophones.  This  is  described  by 
the  formula: 

0(r)  =  £u(rn,|r-rn|/c)|r-rn|,  (1) 

W=*l 

where  n(r„,  t)  means  a  (possibly)  filtered  version  of  u(rn,  t),  that  is, 

«(r „,*)  =  u(Tn,t)*F(t),  (2) 

with  *  denoting  convolution  and  F(t)  a  filter  function.  For  example,  F(t)  might  be 
a  whitening  filter  designed  to  boost  the  higher  frequencies.  In  (1)  the  factor  |r  —  r„| 
is  included  to  compensate  for  signal  loss  due  to  geometric  spreading  of  the  wave 
originating  from  the  point  r.  This  factor  is  optional  and  simulations  show  that  its 
inclusion  does  not  have  a  very  noticeable  effect  on  image  quality.  Note  that  in  (1), 
the  signal  u( r„,  t)  is  evaluated  at  the  time  t  —  |r  —  rn|/c;  this  correctly  describes 
backpropagation  to  the  point  r  if  the  time  origin  (t  —  0)  occurred  when  the  signal 


1 


was  emitted  from  r.  However,  the  time  origin  is  unknown.  This  point  is  discussed 
further  below. 


Now,  for  brevity,  define  u„(r)  =  u( rn,  |r  —  rn|/c)|r  —  rn|,  so  that  (1)  is 


0(r)  =  E“«(r)‘  (3) 

n~l 

The  function  un(r)  is  a  random  variable  with  zero  mean.  Letting  E{-}  denote 
statistical  expectation  (or  ensemble  average),  we  see  that  2£{0(r)}  =  0  since  the 
image  0(r)  is  linear  in  the  field  u(r,  t).  This  precludes  averaging  the  images  on  an 
“amplitude”  basis.  We  therefore  consider  averaging  successive  images  on  an 
intensity  basis,  i.e.,  averaging  |0(r)|2.  This,  however,  leads  to  a  large  “DC  bias”  in 
the  image  that  builds  up  over  time.  An  unbiased  image  estimator  is  obtained  by 
subtracting  the  “DC  component,”  as  follows: 


fl  M  |2j 

f  N 

0(r)  =  £{£*»(*)  [ 

-  ^  EMr)|2  , 

(4) 

(Jn=l  1  J 

ln=l  ) 

which  can  also  be  written 

0(r)  =  J2Y1  (5) 

n=l  m— 1 
n^m 

where  *  means  complex  conjugate.  For  generality,  we  allow  the  fields  ttn(r)  to  be 
complex,  because  we  may  choose  u(r,  f)  to  represent  either  a  real  quantity  or  its 
analytic  signal,  which  is  complex.  Equation  (5)  shows  that,  if  there  is  no  correlation 
between  distinct  traces,  u„(r),  then  the  image  0(r)  vanishes  as  desired.  The  first 
term  in  (4)  bears  some  resemblance  to  the  coherence  measure  called  the  semblance 

W. 


Suppose  we  denote  by  O^(r)  the  fc-th  image  realization.  We  wish  to  average  many 
such  realizations  to  approximate  the  ensemble  average  indicated  by  the  expectation 
operator  £{•}  in  (4).  The  index  fc  may  be  thought  of  as  a  time  index.  At  each  time 
increment,  we  compute  a  complete  image,  which  requires  a  summation  over  all 
geophones  for  each  pixel  in  the  image.  We  can  average  these  images  recursively,  as 
follows.  The  average  of  M  images  is 


i  M 

<3«(r)  = 

1V1  fc=l 

where  the  fc-th  realization  is  defined  by 


0<*>(r)  = 


E  *?>(') 


|n=l 


(2 


N 


E 


(6) 

(7) 


2 


One  can  update  the  average  (6)  recursively  as  more  data  become  available.  That  is, 
given  the  average  formed  after  M  image  realizations,  i.e.,  Ojv/(r),  then  the  updated 
mean  after  the  M  + 1  realization,  0(M+1)(r),  becomes  available  is 


<W*>  =  STT0"(r)  +  MTTo<"+1)(r)'  (8) 

In  practice,  it  may  be  easier  to  record  a  finite  record  of  data,  say,  of  K  samples,  and 
then  reconstruct  K  images  and  simply  average  these,  rather  than  to  use  the  above 
recursive  algorithm  at  every  time  increment.  The  recursive  procedure  could  then  be 
used  on  the  sets  of  averaged  K  images.  If  K  is  not  too  large,  the  latter  approach 
shouldn’t  require  storage  of  an  unmanageable  amount  of  data.  We  can  write  this 
averaging  process  with  respect  to  k  more  explicitly  as  follows.  To  average  K  images 
at  K  times  t*,  we  can  write  the  first  term  in  (4)  as  follows: 


0?\r)  =  ££ 


K 


K 


and  the  second  term  in  (4)  as 


k=  i 


N 


^  (rn , tfc  -f"  |r  rn|/c)  |r  rn 

71—1 


(9) 


0?>(  r)  =  ii;i:Wr„(t  +  |r-r„|/c)|2|r-r„|!.  (10) 

K  fc=ln=l 

Essentially,  we  are  reconstructing  images  (via  coherent  summation  over  all 
geophones  in  (9))  at  all  time  origins,  t*. 


2.  Programming  Considerations 

To  describe  the  algorithm,  let  Np^ejj  be  the  number  of  pixels,  N9eo  the  number  of 
geophones,  and  N#  the  number  of  time  samples  in  the  current  data  record.  Let 
RCNpfeeiajNpeo)  denote  an  array  dimensioned  as  indicated  that  contains  the 
pre-computed  distances  between  each  pixel  and  each  geophone.  Let  the  data  array 
be  given  by  Data(NK,Ngeo)  and  let  Image(Npia;ejs)  denote  an  array  that  stores  the 
image.  The  latter  array  is  initialized  to  zero  before  starting  the  algorithm.  Let  c  be 
the  velocity  of  sound  and  let  dt  denote  the  sampling  interval.  The  algorithm  then 
reads  as  follows. 


3 


1  do  fc=Njr 

2  do 

3  sumi  =  0 

4  sum2  =  0 

5 

6  do 

7  delay  =  R(i,j)/c 

8  ik  =  round  (delay/dt  +  k) 

9  sumi  =  sumi  +  Data  (ik,j) 

10  sum2  =  sum2  +  Data(ifcj)**2 

11  end 

12 

13  image(i)  =  image(i)  +  sumi*  *2  -  sum2 

14  end 

15  end 

Line  7  computes  the  propagation  delay  between  the  jf-th  geophone  and  the  i-th 
pixel.  Line  8  rounds  this  to  an  integer  index  that  is  used  to  select  the  time  index  ik 
in  the  data  array  Data(ifcj).  Note  that  the  incoherent  average  is  performed  with 
respect  to  the  index  k  (time)  and  the  coherent  average  is  performed  with  respect  to 
the  index  j  (geophone). 

3.  Caveats 

This  algorithm  is  very  simple.  It  continually  backpropagates  as  new  data  arrive  and 
sums  the  images  on  an  intensity  basis.  As  described,  the  algorithm  neglects 
variations  in  the  sound  velocity,  assumes  one  type  of  (scalar)  wave,  assumes  perfect 
knowledge  of  the  geophone  coordinates,  and  assumes  that  the  sources  are 
broadband  and  spatially  incoherent  (uncorrelated).  Also,  the  geophone  spacing  is 
assumed  small  enough  so  that  spatial  frequency  aliasing  is  not  a  problem  at  the 
higher  frequencies.  Departures  from  any  of  these  assumptions  will  of  course  tend  to 
degrade  image  quality. 
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Abstract — An  analysis  of  a  passive  seismic  method  for  subsur¬ 
face  imaging  is  presented,  in  which  ambient  seismic  noise  is  em¬ 
ployed  as  the  source  of  illumination  of  underground  scatterers. 
The  imaging  algorithm  can  incorporate  new  data  into  the  image 
in  a  recursive  fashion,  which  causes  image  background  noise  to 
diminish  over  time.  Under  the  assumption  of  spatially-incoherent 
ambient  noise,  an  analytical  expression  for  the  point-spread  func¬ 
tion  of  the  imaging  algorithm  is  derived.  The  point-spread  function 
(PDF)  characterizes  the  resolution  of  the  image,  which  is  a  function 
of  the  receiving  array  length  and  the  ambient  noise  bandwidth.  Re¬ 
sults  of  a  Monte  Carlo  simulation  are  presented  to  illustrate  the 
theory. 

Index  Terms — Acoustic  imaging,  acoustic  tomography,  geophys¬ 
ical  inverse  problems,  geophysical  signal  processing,  geophysical 
tomography,  image  reconstruction,  seismic  inverse  problems. 

I.  Introduction 

N  time-exposure  photography,  low  intensity  light  is  inte¬ 
grated  over  time  to  enhance  image  brightness.  A  similar 
concept  is  applicable  in  acoustics,  wherein  acoustic  or  seismic 
imaging  is  performed  using  ambient  acoustic  noise  as  a  constant 
source  of  “illumination.”  The  idea  of  passive  acoustic  imaging 
was  recently  demonstrated  in  an  ocean  environment  by  Buck¬ 
ingham  and  collaborators,  in  which  submerged  scatterers  were 
imaged  using  ambient  acoustic  noise  generated  by  surface  wave 
action  [  1  ]— [3].  The  authors  refer  to  this  method  as  “acoustic  day¬ 
light  imaging.'*  We  propose  extending  the  concept  to  seismic 
imaging  of  subsurface  scatterers,  including  natural  and  man¬ 
made  underground  structures.  Sources  of  ambient  seismic  noise 
are  microseisms  and  surface  noise  due  to  weather  and  human 
activity. 

Ambient  imaging  is  potentially  attractive  since  it  avoids  the 
high  cost  of  the  deployment  and  operation  of  active  seismic 
sources.  Moreover,  ambient  illumination  is  everpresent,  thus 
permitting  the  “integration**  of  ambient  data  over  an  extended 
time  period.  In  this  way,  relatively  weak  scatterers  should  even¬ 
tually  reveal  themselves  after  a  sufficient  accumulation  of  data. 

The  basis  of  an  imaging  method  using  ambient  illumination 
can  be  explained  as  follows.  Seismic  or  acoustic  waves  im¬ 
pinging  on  a  scatterer  will  create  a  scattered  wavefront  that  will 
manifest  itself  as  a  spatially  correlated  signal  over  a  receiving 
array  on  the  surface.  By  coherently  summing  this  correlated 
signal  using  a  suitable  migration  (or  backpropagation)  proce¬ 
dure,  the  image  of  a  weak  scatterer  should  grow  stronger  over 
time  relative  to  the  uncorrelated  background  noise.  One  can  con¬ 
ceive  of  devising  an  imaging  algorithm  that  develops  an  image 
of  the  target  as  new  data  are  incorporated  into  the  image  in  a 
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recursive  manner.  In  this  way.  storage  of  large  quantities  of  data 
prior  to  forming  the  image  can  be  avoided. 

In  the  seismic  domain,  a  few  previous  attempts  have  been  re¬ 
ported  on  passively  detecting  large  reflecting  structures  (e.g., 
a  bedrock  interface)  [4]-[7].  As  far  as  we  are  aware,  however, 
no  attempts  have  been  made  to  image  scatterers  or  diffracting 
structures  (i.e.,  structures  comparable  to  the  wavelength  of  the 
illuminating  waves)  using  ambient  illumination.  A  more  com¬ 
plete  review  of  published  work  on  passive  seismic  methods  can 
be  found  in  a  recent  report  from  the  Oak  Ridge  National  Lab¬ 
oratory,  Oak  Ridge,  TN,  [8]  that  presents  results  of  a  prelimi¬ 
nary  study  on  the  feasibility  of  subsurface  ambient  imaging.  In 
that  study,  no  imaging  was  attempted  using  real  data,  but  some 
important  facts  were  established.  First,  ambient  noise  in  the  ap¬ 
propriate  frequency  range  (up  to  100  to  200  Hz)  has  adequate 
energy  to  be  used  for  subsurface  imaging,  and  second,  cross-cor¬ 
relation  of  noise  recorded  by  geophones  in  a  two-dimensional 
(2-D)  array  showed  peaks  that  relate  to  seismic  energy  prop¬ 
agating  across  the  array.  It  was  determined  that  the  source  of 
this  spatially-coherent  energy  was  traffic  from  a  nearby  road. 
An  implication  of  the  report  is  that  the  ability  to  suppress  sur¬ 
face  wave  energy  by  digital  filtering  or  optimal  array  design  is 
critically  important  in  separating  ground  roll  from  the  upward 
propagating  body  waves  scattered  from  subsurface  structures. 

An  essential  difference  between  active  and  ambient  illumina¬ 
tion  is  that,  in  the  latter  case,  no  time  origin  exists  with  respect 
to  which  propagation  delays  can  be  referenced.  The  illumina¬ 
tion  is  effectively  continuous  and  random.  A  greatly  simplified 
analysis  is  possible  if  the  illuminating  waves  can  be  regarded  as 
spatially  incoherent.  Spatial  incoherence  is  an  idealization  im¬ 
plying  that  the  wave  fields  scattered  from  two  distinct  points  in 
space  are  statistically  uncorrelated  over  time.  If  waves  are  ar¬ 
riving  from  one  direction  from  a  distant  and  localized  source, 
the  assumption  of  spatially-incoherent  illumination  will  break 
down.  However,  if  the  source  of  noise  is  primarily  from  the 
surface,  which  is  generally  the  case  at  higher  frequencies,  then 
the  assumption  of  spatial  incoherence  is  probably  quite  good 
when  the  averaging  is  performed  over  an  extended  time  period. 

It  can  be  shown  that  the  spatial  resolution  in  an  image  derived 
from  ambient  noise  is  determined  essentially  by  two  factors:  the 
length  of  the  geophone  array  and  the  spatial  coherence  length  of 
the  noise.  The  latter  is  approximately  the  speed  of  sound  divided 
by  the  noise  bandwidth.  For  analytical  tractability,  we  assume 
scatterers  embedded  in  a  background  of  constant  velocity.  Al¬ 
though  the  uniform  background  assumption  is  somewhat  artifi¬ 
cial,  our  primary  aim  in  this  paper  is  to  investigate  the  resolution 
limits  of  a  passive  imaging  scheme  and  relate  it  to  the  statistical 
properties  of  the  seismic  noise. 

Below,  an  expression  for  the  point-spread  function  (PSF, 
which  is  the  image  of  a  point  scatterer)  of  an  ambient  imaging 
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system  is  derived  assuming  complete  spatial  incoherence  of  the 
noise.  Suppose  that  p{v\ :  t\ )  and  p( r2 :  !•>)  are  signals  scattered 
from  the  point  rj,  and  r>-  Then,  complete  spatial  incoherence 
can  be  expressed  as 

E[p(rx:tx)p{\'2\hY]  =  -  M'Hn  “  L>)  0) 

where  £[•]  denotes  statistical  expectation  (i.e..  an  ensemble  av¬ 
erage),  <5( r i  -  r2)  is  a  three-dimensional  (3:D)  Dirac  delta  func¬ 
tion.  and  P(t)  is  the  temporal  autocorrelation  function  of  the 
noise.  Here,  we  assume  a  stationary  random  process,  implying 
that  the  argument  of  P(t)  is  a  function  of  the  time  difference 
T  _  tx  _  t2.  In  ( I),  *  denotes  complex  conjugate,  since  we  re¬ 
serve  the  option  of  working  with  analytic  signals.  We  should 
point  out  that  the  relation  (1)  makes  sense  for  a  spatially  inco¬ 
herent  source,  but  for  a  scattering  problem.  (I)  will  always  rep¬ 
resent  an  approximation  whenever  the  autocorrelation  function 
P(t)  has  some  finite  duration  (i.e.,  when  the  noise  is  bandlim- 
ited).  This  is  because  an  incident  waveform  that  is  correlated 
over  a  finite  time  (given  by  the  reciprocal  of  its  bandwidth)  will 
give  rise  to  scattering  that  is  spatially  correlated  over  distances 
up  to  the  spatial  correlation  length  of  the  noise,  given  by  c/B, 
where  D  is  the  bandwidth  of  the  noise,  and  c  is  the  speed  of 
sound.  As  a  consequence,  the  Dirac  delta  function  in  ( 1 )  will  be 
broadened  by  this  amount.  If  we  consider  scatterers  farther  apart 
than  c/B.  then  the  relation  (I)  effectively  holds.  In  an  imaging 
algorithm,  it  thus  makes  sense  to  make  our  pixel  size  larger  than 
this  correlation  length.  Image  resolution  can  also  be  improved 
to  some  extent  by  applying  a  prewhitening  filter  to  the  received 
signals  prior  to.  image  formation.  Such  a  filter  can  help  flatten 
the  spectral  response  over  the  available  bandwidth  which,  in  the 
spatial  domain,  translates  into  a  narrow  point  response. 

The  assumption  of  complete  incoherence,  expressed  by  (1), 
considerably  simplifies  the  derivation  of  an  analytical  expres¬ 
sion  for  the  FSF.  Later,  we  compare  the  analytically  derived  PSF 
to  a  more  realistic  computation  of  the  PSF  based  on  a  Monte 
Carlo  simulation  of  bandlimited  stochastic  noise  using  a  random 
number  generator.  We  shall  see  that  the  Monte  Carlo  simulation 
shows  good  agreement  with  the  analytically  derived  PSF.  The 
simulated  images  of  several  scatterers  were  obtained  for  var¬ 
ious  numbers  of  exposures,  where  each  “exposure  refers  to  one 
image  realization  generated  from  one  set  of  independent  data. 
“Independent  data”  means  that  the  interval  of  time  separating 
data  sets  is  taken  to  be  greater  than  the  reciprocal  of  the  band¬ 
width  of  the  noise.  Assuming  a  bandwidth  extending  to  200  Hz, 
this  implies  an  exposure  interval  of  0.005  s.  Thus,  1 00  000  expo¬ 
sures  would  take  no  shorter  than  about  500  seconds  of  recording 
time.  We  will  see  that,  as  the  number  of  exposures  that  are  av¬ 
eraged  increases,  a  gradual  reduction  in  the  uncorrelated  back¬ 
ground  noise  in  the  image  is  evident  as  data  are  continuously 
incorporated. 

In  the  analysis  that  follows,  we  make  several  idealizing 
assumptions  for  the  sake  of  analytical  tractability,  such  as 
the  assumption  that  stochastic  scattering  can  be  modeled 
as  a  weighted  superposition  of  incoherent  sources,  and  that 
multiple  scattering  is  sufficiently  small  to  neglect.  It  is  likely, 
however,  that  the  success  of  the  method  will  not  rest  crucially 
on  any  of  these  idealizations,  although  performance  may  depart 
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Fig.  I.  Point-spread  function  (PSF)  for  a  scatterer  at  a  depth  of  30  m.  computed 
using  ( 16)  with  surface  data  only. 

somewhat  from  the  results  predicted  by  the  analysis  based  on 
our  simplifying  assumptions. 

II.  Imaging  Algorithm 

We  model  the  earth  as  a  distribution  of  scatterers  of  varying 
strengths  (scattering  cross-sections)  embedded  in  a  uniform 
background  of  constant  velocity  r.  The  extension  to  the  case 
of  nonuniform  velocity  can  be  made,  provided  that  velocity 
information  can  be  derived  by  independent  means.  We  assume 
that  the  scattering  medium  is  illuminated  with  broadband, 
spatially-incoherent  noise.  The  noise  is  then  passively  recorded 
on  the  surface  of  the  earth  at  N  measurement  points  over 
some  interval  of  time.  Our  objective  is  to  reconstruct,  or 
image,  the  scattering  density  distribution  by  backpropagating 
(migrating)  the  recorded  noise  into  the  earth.  Because  the 
source  of  illumination  is  always  “turned  on,”  an  image  of  the 
scattering  distribution  can  be  gradually  improved  over  time 
by  recursively  incorporating  new  data.  The  imaging  algorithm 
coherently  sums  over  the  wavefront  scattered  from  a  point  in 
the  earth  (equivalent  to  coherently  summing  along  a  hyperbolic 
path  in  the  space-time  domain)  and  assigns  this  value  to  the 
corresponding  image  point.  This  process  is  repeated  for  all 
image  points  until  an  entire  image  is  generated.  Multiple 
images  created  over  successive  time  intervals  are  then  added 
on  an  intensity  basis  as  described  below.  As  more  data  are 
incorporated,  simulations  show  that  the  stochastic  variability  in 
the  (mean)  image  due  to  the  random  nature  of  the  illumination 
gradually  diminishes,  and  weak  scattering  features  become 
more  apparent.  After  a  certain  large  number  of  iterations,  the 
image  ceases  to  improve  since  all  image  noise  has  effectively 
been  averaged  out.  Any  remaining  artifacts  (e.g..  the  finite  PSF 
width  and  sidelobe  structure)  are  due  to  the  inherent  resolution 
limits  of  the  system  imposed  by  finite  aperture  and  bandwidth 
constraints. 
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Fig.  2.  Normalized  images  of  three  scattered  using  stochastic  data  for  different  numbers  of  exposures  with  surface  data  only. 


Below,  we  describe  the  imaging  algorithm  and  derive  from  it 
the  PSF  of  the  system,  which  may  be  defined  as  the  image  of 
a  single  point  scatterer.  The  PSF  characterizes  the  lateral  and 
depth  resolution  of  the  imaging  system. 

Let  u(rn,t)  denote  the  field  recorded  at  the  point  r,t  on  the 
Earth's  surface.  Then  the  image  generated  by  the  process  of 
backpropagation  is  given  by 


0{ r)  =  4tt  «(«•>■■  |r  -  r, ,(/<:)  |r  -  r„| 

M  =  1 


where  u(rn,  t)  is  a  filtered  version  of  u(rtl1 1).  That  is 

u(rnJ.)  =  u(rnjt)  *  F(t)  (3) 

with  *  denoting  convolution  and  F(t)  a  specified  filter  function. 
For  example,  F(t)  could  be  a  whitening  filter.  In  (2),  the  factor 
47r|r  -  rn  |  is  included  to  compensate  for  signal  loss  due  to  geo¬ 
metric  spreading  of  the  wave  originating  from  the  point  r. 

For  brevity,  define  un( r)  =  47ru(ra,  |r  -  rn|/c)  |r  -  rfl|  so 
that  (2)  may  be  written 

.v 

O(r)  =  'jr  nj r). 

n=  I 


(2) 


(4) 
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The  noise  field  «„(r)  is  a  random  variable  with  zero  mean.  Let¬ 
ting  £{•}  denote  statistical  expectation  (or  ensemble  average), 
we  see  that  E{0( r)}  =  0.  since  the  image  0( r)  is  linear  in  the 
field  ii(r.  (.).  This  precludes  averaging  the  images  on  an  “ampli¬ 
tude”  basis.  We  therefore  consider  averaging  successive  images 
on  an  intensity  basis  (i.e.,  averaging  |C?(r)|2).  This,  however, 
leads  to  a  large  DC  bias  in  the  image  that  builds  up  over  time. 
An  unbiased  image  estimator  is  obtained  by  subtracting  the  DC 
component  as  follows: 


O(r)  =  E 


E 


jV 

E  Kt»(r)l2 


(5) 


which  can  also  be  written 

N  N 

=  (6) 

n=l  m  =  l 
rn 


Equation  (6)  shows  that  if  there  is^no  correlation  between  dis¬ 
tinct  traces  «„(r),  then  the  image  O(r)  vanishes  as  desired.  Oth¬ 
erwise,  it  is  positive.  The  first  term  in  (5)  bears  some  resem¬ 
blance  to  the  coherence  measure  called  the  semblance  [9]. 

Suppose  we  denote  by  0{k)( r)  the  A;th  image  realization.  We 
wish  to  average  many  such  realizations  to  approximate  the  en¬ 
semble  average  indicated  by  the  expectation  operator  £{*}  in 
(5).  That  is,  the  average  of  M  images  is 


0.u(r )=T  E  <?<fc)(r 


(7) 


k=i 


where  the  A:th  realization  is  defined  by 


0(k\ r) 


iV 


n- 1 


E  i^wi2.  (8) 


n=l 


To  avoid  storing  large  amounts  of  data,  one  can  update  the  av¬ 
erage  (7)  recursively  as  more  data  become  available.  Given  the 
average  Om{ r)  formed  after  M  image  realizations,  the  updated 
average  (9ty/+l(r)  after  the  M  -f  1  realization  0^AI+l^{r)  be¬ 
comes  available  is 
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Fig.  3.  PSF  for  a  scniterer  computed  using  { 16).  assuming  both  surface  and 
borehole  data. 


is  to  reconstruct  a(r).  Suppose  random  noise  is  illuminating  the 
medium  from  many  directions,  and  let  p( r;  t)  denote  the  noise 
scattered  from  the  point  r.  Then  the  signal  received  at  a  point 
rtt  on  the  array  may  be  written  as  the  following  integral  over  the 
medium 


"■(r  ,..()=  IJj 


p[r':/  -  jr'  -  r„|/c) 
47r|r'  -  r„| 


r/V.  (10) 


In  this  simple  model,  multiple  scattering  has  been  neglected.  For 
a  3-D  reconstruction,  we  record  data  over  a  2-D  array  on  the 
surface  of  the  earth.  Letting  r(l  represent  a  point  on  this  array 
and  .4  the  region  of  the  surface  occupied  by  the  array,  then  the 
backpropagation  operation  can  be  expressed  as  follows,  where, 
for  generality,  we  again  allow  filtering  prior  to  backpropagation: 


0{r)  =  4tt  Jj  d7ra  u(ra.  |r  -  rfi|/c)  |r  -  r„|  (11) 

.4 

which  is  the  continuous  space  analogue  of  (2).  Once  again, 
u(rnJ.)  is  a  filtered  version  of  the  recorded  data  u(r„,  t).  i.e.. 


<?A/+l(r)  = 


M  O.u(r)  + 


1 


M  +  1 


M  +  1 


C?(-u+1>(r) 


(9) 


III.  Derivation  of  System  Point  Spread  Function 


n(rnJ.)  =  u(ruJ.)*F(L)  (12) 

for  a  specified  filter  function  F{t). 

To  compute  the  image,  we  substitute  ( 1 0)  into  (11)  and  inter¬ 
change  orders  of  integration.  This  gives 


Below,  we  derive  an  analytical  form  for  the  PSF.  Our  basic 
objective  is  to  evaluate  analytically  the  first  term  in  (5),  which 
represents  the  actual  image.  The  purpose  of  the  second  term 
is  essentially  to  eliminate  a  DC  bias,  or  constant  background, 
from  this  image.  For  analytical  tractability,  we  assume  that  the 
receiving  array  is  continuously  sampled,  in  which  case,  the  sum 
over  n  in  (5)  becomes  an  integral  over  the  array  aperture.  We 
model  the  medium  as  a  continuum  of  scatterers  embedded  in  a 
background  of  constant  velocity  c.  The  density  of  scatterers  at 
the  point  r  is  denoted  by  rr(r):.  That  is,  <r(r)  can  be  regarded  as  a 
continuously  varying  scattering  cross-section,  and  our  objective 


0(r)  =  ffl  i/V< 7(r')  JJ  d2 r„ 

A 

■  ^  p[r'\  |r  -  r„|/<:  -  |r'  -  ru|/c]  (13) 

lr  -  r«l 

where  p{rj.)  =  p(v:t)  *  F(t ).  Now  take  the  absolute  square 
of  (13),  perform  an  ensemble  average,  and  use  the  relation  (1). 
This  gives  rise  to  the  following  result: 

£(|(?(r)|2]  =  IJJ  d3r'  cj(r')2  PSF(r:r')  (14) 
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(a)  1  exposure 


(b)  10  exposures 
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Fig.  4.  Normalized  images  of  three  scatterers  using  stochastic  data  for  different  numbers  of  exposures  assuming  both  surface  and  borehole  data. 
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where  F(w)  is  the  Fourier  transform  of  the  filter  function  F(t.). 
Substituting  this  into  (15)  results  in 

PSF(r;  r')  =  f  duj  P(u>)  |F(w)|2|/L,(r;  r')|2  (16) 

J  —  OC 

where 

■4“<r;r')s/ 

•  exp['iu/(|r  -  r«|  -  |r'  -  ra\)/c\.  (17) 

Clearly,  at  each  frequency  w,  the  function  .4W( r;  r')  will  become 
large  when  the  image  point  r  coincides  with  the  scattering  point 
r'  since  then  the  exponent  vanishes.  In  the  next  section,  we  show 
a  numerical  computation  of  the  PSF  using  (16)  and  (17),  and 
compare  it  to  a  Monte  Carlo  simulation. 

IV.  Simulations 

Fig.  1  shows  an  image  of  PSF  PSF(r;  ro)  obtained  by  nu¬ 
merically  evaluating  (16)  and  (17)  over  a  linear  array.  In  com¬ 
puting  the  PSF,  the  scattering  point  was  placed  at  r'  =  (0.  -30) 
in  units  of  meters.  Twenty  receiving  elements  were  assumed, 
spaced  at  5-m  intervals  (implying  an  array  length  of  95  m).  The 
power  spectrum  P{u)  was  assumed,  tor  simplicity,  to  extend 
uniformly  to  200  Hz.  A  velocity  of  sound  of  500  m/s  was  used. 
We  note  that  if  the  velocity  is  scaled  up  and  the  noise  bandwidth 
remains  unchanged  then  the  appearance  of  all  images  will  re¬ 
main  unchanged,  provided  all  distances  are  scaled  in  proportion 
to  the  velocity. 

Fig.  2  shows  the  results  of  a  Monte  Carlo  simulation  using  (7) 
and  (8),  which  combine  multiple  image  realizations,  each  gener¬ 
ated  by  the  backpropagation  process  defined  by  (2).  In  the  sim¬ 
ulation,  three  scatterers  were  assumed  at  coordinates  (-12.5, 
-20),  (-2.5,  -35),  and  (12.5,  -45).  We  regarded  each  of  these 
three  points  as  behaving  effectively  as  independent  sources  of 
broadband  random  noise.  As  noted,  this  is  a  reasonable  model 
when  the  noise  correlation  distance  is  smaller  than  the  scatterer 
separations.  Finally,  surface  interactions  and  multiple  scattering 
between  the  scatterers  were  neglected  in  the  simulation. 

The  number  of  exposures  indicated  in  Fig.  2  is  defined  as  the 
number  of  terms  M  in  the  summation  (7).  In  each  image  real¬ 
ization,  backpropagation  to  a  10  x  10  array  of  pixels  was  per¬ 
formed  with  a  center-to-center  pixel  separation  of  5  m.  The  sim¬ 
ulated  scattered  waveform  ti(r,  t)  was  created  using  a  random 
number  generator  to  produce  a  zero-mean,  uniformly  distributed 
sequence  of  samples  at  a  sampling  interval  of  0.0025  s,  corre¬ 
sponding  to  a  Nyquist  rate  of  200  Hz.  This  implies  a  noise  corre¬ 
lation  length  c/ B  equal  to  2.5  m  when  B  =  200  Hz  and  c  =  500 
m/s.  Thus,  the  center-to-center  pixel  spacing  is  greater  than  the 
correlation  length.  The  numbers  of  exposures  averaged  are  indi¬ 
cated  in  the  figure.  Each  image  has  been  normalized  to  unity  to 
show  clearly  how  background  artifacts  are  reduced  as  the  expo¬ 
sure  number  increases.  In  this  particular  simulation,  the  image 


quality  reaches  a  steady  state  above  about  1000  exposures.  This 
means  that  nearly  all  of  the  stochastic  variability  has  been  elim¬ 
inated  at  this  point,  and  the  remaining  artifacts  (PSF  broadening 
and  sidelobe  structure)  are  determined  by  the  inherent  band¬ 
width  and  aperture  limitations  of  the  system. 

The  above  simulation  assumed  an  array  of  20  geophones  on 
the  surface  of  the  earth.  Another  simulation  was  performed  with 
both  surface  geophones  and  borehole  geophones.  The  same 
scatter ers,  velocity,  and  noise  bandwidth  were  assumed.  In  this 
case,  however,  60  geophones  were  assumed  in  the  simulation, 
20  on  the  surface,  20  distributed  vertically  down  a  borehole  to 
the  left  of  the  scattering  region,  and  20  in  a  borehole  to  the  right 
of  the  scattering  region.  The  three  geophone  array  segments, 
one  surface  and  two  borehole,  were  assumed  to  be  of  the  same 
length  (95  m).  The  PSF  for  this  array  configuration  was  again 
computed  using  (16)  and  (17)  and  is  shown  in  Fig.  3.  Next,  the 
results  of  a  Monte  Carlo  simulation  are  shown  in  Fig.  4.  Note 
that  both  the  PSF  and  the  simulated  images  are  sharper  and 
more  symmetrical  due  to  the  extended  system  aperture  created 
by  the  addition  of  the  boreholes  data. 


V.  Conclusion 

A  theoretical  investigation  of  the  possibility  of  imaging  un¬ 
derground  scatterers  employing  ambient  illumination  was  pre¬ 
sented.  For  simplicity,  we  assumed  that  the  scatterers  were  em¬ 
bedded  in  a  constant  velocity  background,  but  the  extension  to 
nonuniform  velocity  is  expected  to  yield  qualitatively  similar 
results.  Under  the  assumption  of  spatially  incoherent  noise,  an 
analytical  expression  for  the  PSF  was  derived.  This  expression 
characterizes  both  transverse  and  depth  resolution  and  defines 
their  dependence  on  array  length  and  the  ambient  noise  band¬ 
width.  Monte  Carlo  simulations  were  performed  to  demonstrate 
that  multiple  exposures,  added  together  as  data  are  continuously 
acquired,  will  diminish  uncorrelated  background  noise  and  im¬ 
prove  image  quality. 
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